package scu.maqiang.matlabcases;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T61;
import scu.maqiang.mesh.Mesh2T6;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class PoissonOnUnitDisk {

	public static void main(String[] args) {
		Mesh2T6 mesh2D = new Mesh2T6();
		mesh2D.readMesh(MeshFormat.msh, "matlabMeshes\\finalnitdisk.msh");
		FES2T61 fs = new FES2T61(mesh2D);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		double[] x = new double[fs.GetNdof()];
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.assembleSource(new double[] {1.0}, BVPType.COMMON, RHS);
		fs.applyBC_MBN(A, RHS, Direct.All, 0.0, 1, 2, 3, 4);
		NewIterSSolver solver = new NewIterSSolver(A);
		solver.PCGSSOR(RHS, x, 1.65, 1);
		ScalarFunc exactFunc = (xy, label, param) -> 0.25 * (1 - xy[0] * xy[0] - xy[1] * xy[1]);
		double[] exactSolution = fs.valueFromFunc(exactFunc, null);
		double[] error = MVO.add(x, -1.0, exactSolution);
		System.out.println("Error min: " + MVO.min(error) + "\tmax: " + MVO.max(error));
	}
}
